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I review recent progress on algorithms for calculating quark propagators and for simulating full QCD. 
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For the sake of brevity and consistency, I ex- 
cluded from this review many interesting papers 
which discussed algorithms not directly relevant 
for QCD. I apologize to their authors, and refer 
the reader to their contributions in this volume. 

The bottleneck of quenched lattice QCD is the 
calculation of quark propagators; for full QCD 
simulations, the Monte Carlo algorithm itself is 
the bottleneck. I discuss these two issues in Sec- 
tions 1 and 2. 

1. Quark propagator calculation 

The task is to solve the linear system 
A{U,mg)x = b, where A{U,mg) is the discretized 
Dirac operator. This task is to be repeated many 
times, over matrices A constructed from different 
gauge fields U with statisticahy similar proper- 
tics, usually for a whole set of right-hand sides 
b, and often for a whole range of quark masses 
rrig. Let us first consider the solution of a single 
system. 

The matrix A is not hermitian, so that the 
standard Conjugate Gradient (CG) algorithm can 
only be applied to the system A^Ax = A^b (it 
is then abbreviated CGNE: CG on the Normal 
Equations). CGNE is expected to have mediocre 
convergence properties, because the condition 
number of A^A is the square of that of A. There- 
fore, until recently, the preferred method was 
nearly the simplest: do a line minimization of 
II Ax — 6 II at each iteration, in the search di- 
rection defined by A times the previous resid- 
ual r = Ax — b, but orthogonal to the previous 
search direction. This algorithm, called Conju- 
gate Residual (CR) [also cahed CR(1) or GM- 



RES(l)] offers no minimum rate of convergence, 
and is therefore supplemented by CGNE: one 
switches to CGNE when convergence of CR be- 
comes unsatisfactory. 

The situation has changed with the recognition 
that there exist other iterative solvers, designed 
for non-hermitian matrices, which converge faster 
than CR. Early work on this subject [1,3] stressed 
the advantages of BiCG-type methods, and quan- 
tified the gain over CR and CG: the number of 
matrix- vector products is reduced by a factor ~ 2. 
The question then arises: can one do better ? 

It is possible to answer this question for all cur- 
rent algorithms because they are all Krylov meth- 
ods. By that one means that at each iteration k, 
an approximate solution x^ is found in the vector 
space (Krylov space) 

£k = span{ro, Aro, A'^ro, A'^-'^ro} (1) 

Different algorithms come more or less close to 
finding in the vector x^ which minimizes the 
norm of the residual || Ax — b \\. An exact de- 
termination of Xk requires the construction of an 
orthonormal basis of £k- When A is hermitian 
positive definite, CG has the remarkable prop- 
erty of building this basis recursively, with only 
2 inner products per iteration, and finding Xk ex- 
actly: it can be called an "optimal" algorithm, 
in the sense that the number k of matrix-vector 
products required to reduce the norm of the resid- 
ual to a given tolerance is the minimum possible. 
When A is non-hermitian, the progressive con- 
struction of an orthonormal basis of £k for finding 
Xk necessitates the evaluation of k inner products 
< ej|ej > at each iteration k, and the storage of k 
vectors: the work grows like fc^ and the memory 
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Figure 1: Typical eigenvalue spectra for the Dirac 
and for staggered (KS) fermions 

requirement like k. These two demands render 
this algorithm, called (full) GMRES (Generalized 
Minimum RESidual) impractical except for tests. 

But GMRES will achieve for non-hermitian 
matrices A what CG did in the hermitian case: 
it will provide a lower bound on the number 
of matrix-vector products necessary to reach the 
stopping criterion. We can then assess how far 
from optimal other algorithms are, by compari- 
son with this lower bound. We do this in turn 
for Wilson and Staggered fermions, then review 
directions for further progress. 

1.1. Wilson fermions 

The matrix is A = 1 — kM, with symmetries: 

• EMS = — M, where S = ±1 on even/odd lat- 
tice sites. This is because M is a hopping matrix, 
connecting nearest neighbours on the lattice {M 
has a block-block-block-tridiagonal structure) . 

• 7bM75 = Aft. This is inherited from the con- 
tinuum operator. 

These symmetries respectively imply that eigen- 
values of M come in opposite pairs, and are real 
or come in complex conjugate pairs (real eigen- 
values only appear in the confined phase) . A typ- 
ical spectrum of A/ 2 is shown in Fig. la (confined 
phase); as /3 decreases to 0, fluctuations in A in- 
crease and the spectral density of M becomes uni- 
form in the disk of radius 4 centered at the origin 
(see Fig. lb, where the spectrum of A/2 is shown). 



atrix A/2, for Wilson fermions at /3 = 5 and /3 = 0, 

A comparison of algorithms at /3 = (Fig. 2) 
shows the advantages of non-hermitian solvers, 
in particular those based on the Bi-Conjugate 
Gradient (BiCG, [6]). BiCG constructs 2 Krylov 
spaces, Sk (l) and 

£k = span{so,AK^o,A^^so, A^''-' sq} (2) 

In these 2 spaces, sequences {ui} and {vj} 
which are mutually orthogonal are built: 
< Uk\vi >= y I < k. This orthogonality prop- 
erty ensures convergence after n steps, n being 
the rank of the matrix, in exact arithmetic, al- 
though convergence need not be monotonic. But 
it also contains the germ of numerical instabil- 
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Figure 2. Norm of the residual versus number 
of matrix-vector products, for various solvers at 
/3 = 0,K = .2475 (4^ lattice). 
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ity: the vanishing of < Uk\vi > is the result of 
many cancellations; and Xk, which is built from 
a linear combination of Auk,Uk, and Uk-i, can 
have large numerical round-off errors. The situa- 
tions of breakdown (division by zero - see below) 
or near-breakdown may occur, and "look-ahead" 
or "stabiHzed" BiCG algorithms are designed to 
handle these difficulties gracefully. 

However, one can try the simplest algorithm, 
which uses the 75-symmetry of M to obtain the 
sequence {vj} for free, by just taking Sq = ^^Tq. 
In pseudo-code, this algorithm, which we call 
BiCG75, reads 

Choose starting solution xq. 
rQ = h — A* Xq 
Po = ro 

So =< ro|7B|ro > 

for n = 0, 1, ... until convergence 

= Sn/ < Pn\l5A\pn > 
X„+i = Xn+UJn*Pn 
rn+1 = rn - UJn* A* Pn 

S„+i =< r„+i|75|r„+i > 

Pn+1 = r„+l + 6„+i /Sn*Pn 

end for 

It is identical to CG, except that a 75 has 
been inserted in the inner products. The pos- 
sibility of breakdown comes in the division by 
< Pn\lhA\pn > and by We have experi- 

enced no such breakdown in solving several 10^ 
systems, on a 64-bit Cray computer, with a max- 
imum lattice size 8'^ x 16; but the reader who uses 
a 32-bit machine would be well-advised to use a 
stabilized variant [4]. As Fig. 2 shows, BiCG7B 
compares equally with BiCGStab2, our best sta- 
bilized BiCG algorithm, but requires fewer stor- 
age vectors and fewer inner products. We like it 
mostly for its simplicity. 

Comparing BiCG75 to GMRES as a function 
of K in Fig. 3 reveals that the number of matrix- 
vector products required by BiCG75 is within 
~ 10% of the absolute minimum, over a compre- 
hensive range of quark masses.^ 
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Figure 3. Number of matrix- vector products 
needed to reduce the norm of the residual by 10^°, 
for various solvers, as k is varied. The quark 



mass TUq is estimated as (k 



^)/2, with 



Ke = .1694. The lattice size is 8^ x 16; /3 = 5.7. 

In my opinion, it is not worth looking for a 
better algorithm: BiCG, in its stabilized or in its 
75 version, is quasi-optimal. 

1.2. Staggered fermions 

In this case, the matrix is A = ml +iB, where 

• TjBT, = —B as for Wilson fermions; 

• B^ = B. 

Thus the eigenvalues of B are real and come in 
opposite pairs. The spectrum of A lies on a ver- 
tical line segment (Fig.lc). 

We compared GMRES with CGNE. They re- 
quire the same number of matrix-vector prod- 
ucts, exactly proportional to 1/m. Thus CGNE 
is optimal, confirming the common wisdom that 
nothing can be done to improve convergence over 
CGNE. The question is why. 

A theorem by Voevodin [5] states that the or- 
thogonahzation of the basis in GMRES, which 
normally requires 0{k) inner products at the A;*'' 
iteration, can be reduced to a short recurrence 
and 0(1) inner products whenever matrix A is of 
the form 



A 



\B + al), 61 e 7^, creC 



(3) 



^Admittedly, the comparison sliould be repeated on an 
ensemble of gauge configurations larger than 8^ x 16. The 
cost of GMRES has so far prevented us from doing so. 



Therefore a CG-like algorithm must exist for stag- 
gered fermions. Why is it exactly CGNE ? 
All Krylov methods effectively build at itera- 
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tion k a polynomial Pk{A) in the matrix A, such 
that Pk{A) « A~^, and apply it to the right- 
hand side b. The polynomial built by GMRES, 
Pk{A) or equivalcntly Qk{B), is optimal; and be- 
cause the spectrum of B is even, the polynomial 
Qk{B) will also be even. CGNE builds a poly- 
nomial in A'^ A = m?l + B"^ , i.e. it is an even 
polynomial in B by construction; and it is the 
optimal polynomial in B"^ because CG is optimal 
for a hermitian positive system. Therefore the 
2 polynomials built by GMRES and CGNE are 
both even in B and optimal: they are identical. 

Thus optimal and near-optimal solvers are 
available for staggered and Wilson fermions re- 
spectively. The next improvement should come 
from prcconditioncrs: but it should be clear that 
polynomial-type preconditioners will not reduce 
the amount of work needed, since they do not ex- 
tend the search for a solution outside the original 
Krylov space, which is already explored almost 
perfectly. 

1.3. Further progress 

Further gains can be achieved by avoiding re- 
dundant work during the repetitive solution of 
related systems. Repetition over masses, right- 
hand sides, and gauge configurations are dis- 
cussed in this order. 

1.3.1. Multiple masses 

The structure of matrix A lends itself to an eco- 
nomical calculation of several systems 
(m,l -I- iB)xi = b, i = 1, m 
(using staggered fermions to be specific) , because 
the Krylov space is invariant under a change 
of the quark mass m,. It is sufScient to apply 
the Lanczos process to matrix B once, and then 
to recombine the Lanczos vectors for any set of 
masses. Only one matrix-vector product per it- 
eration is necessary. The savings are machine- 
dependent, but can be considerable, well beyond 
those obtained by the usual procedure of taking 
the last solution as an initial guess for the next 
lighter quark mass. The price to pay is the stor- 
age of all Lanczos vectors of B on disk (for an 
a posteriori reconstruction of other propagators 
for any mass) , or of an extra pair of Lanczos vec- 
tors per mass value in memory (for on-the-fly re- 



construction). Even-odd preconditioning can be 
preserved at an overall cost of 2 in the above. 

Details have been well explained in [2]. This 
simple idea should greatly improve the status of 
quark mass interpolations in quenched QCD. 

1.3.2. Multiple right-hand sides 

Hadron propagators are built from quark prop- 
agators with sources (right-hand sides) spanning 
the color and Dirac spaces. Sources having dif- 
ferent x-space profiles are often desired for varia- 
tional purposes. Thus one is lead to solving 
Axi = bi, i = 1, ...,m 
with m ~ 0{12) or more. 




Figure 4. Logio of the norm of the residual, as a 
function of the number of matrix- vector products, 
for multiple right-hand sides, showing how the 
system "learns" to converge faster. The gauge 
configuration is the same as in Fig. 3; k = .167. 

This is the subject of intense work in the nu- 
merical analysis community [7]. Let me report 
some preliminary results [8]. The idea is to keep 
in memory the most significant search directions 
as a rectangular matrix Q, and then to "deflate" 
the system by solving {l — QQ^)Ax = {l — QQ^)b. 
A nested process allows the system to progres- 
sively refine its set of vectors Q, and to "learn" 
as it iterates. Only one right-hand side is con- 
sidered at a time. Fig. 4 shows the history of the 
norm of the residual as a function of the number n 
of matrix- vector products, for 3 sources on differ- 
ent lattice sites. On the second r.h.s., the system 
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has learned rather well the significant directions 
which slow down convergence, and progresses a 
little more on the third system. Nonetheless, the 
reduction in n is only 0(2). It will increase as the 
right-hand sides become more linearly dependent. 

1.3.3. Multiple configurations 

The eigenvalue spectrum of matrix A, over- 
all, varies little from one configuration {[/} to 
another at the same /3. In fact, it varies little also 
as one scales the size of the lattice: the average 
density of eigenvalues in some region of the com- 
plex plane essentially scales like the volume of the 
lattice. Thus valuable information on the eigen- 
value spectrum, and possibly the structure of the 
eigenvectors after gauge-fixing, could be obtained 
from a small lattice study, e.g. for the purpose of 
building a preconditioner. 

2. Full QCD Monte Carlo 

2.1. Extrapolation methods 

Rather than facing the cost of a full QCD simu- 
lation, one may consider cheaper approximations 
which extrapolate to the full QCD result. The 
simplest extrapolation is in n/ , the number of fla- 
vors. This approach, considered a long time ago 
[9] , has been revived under the name of "bcrmion- 
s" [10]. Consider adding to the gauge action a 
term (f)^l — KM)^1 — K,M)(f), ie. a bosonic term 
with fcrmion-likc interaction (hence the name). 
Integration over gives det^'^(l — kM), as would 
result from -2 quark flavors. Additional pairs of 
bermions can, if desired, mimic n/ = —4,-6,.... 
Together with the quenched result n/ = 0, one 
can extrapolate to n/ > 0. A linear extrapolation 
appears surprisingly good for heavy quarks, when 
compared with full ri/ = 2 results. Bermions are 
attractive because they are cheap to simulate; the 
systematic error in the n/ extrapolation is un- 
clear. 

Another strategy proposed in [12] consists in 
applying to quenched results an estimate of the 
correction due to unquenching. The leading eflPect 
of dynamical fermions is a renormalization of the 
gauge coupling, which can be incorporated in the 
quenched Monte Carlo at no cost. The next-order 
correction is then computed for each observable. 
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Figure 5. Static potential in the Schwinger model. 
The white circles are obtained from quenched re- 
sults corrected for unquenching. From [11]. 

The results are remarkable: the remaining differ- 
ence with full dynamical quarks is buried in the 
statistical noise. This is illustrated in Fig. 5, taken 
from [11], for the static potential in the Schwinger 
model. This is yet another conflrmation that dy- 
namical fermions generate an essentially local ef- 
fective action, unless they are very Hght. It makes 
one wonder anew about possibilities of guiding 
a full QCD simulation with a short-range effec- 
tive action, and enforcing correct sampling with a 
Metropolis test: the acceptance will fall exponen- 
tially with the volume but, on the other hand, the 
autocorrelation time will be reduced by a large 
factor, and probably also a power of the volume, 
over Hybrid Monte Carlo. For moderate volumes 
and quark masses, this simple strategy may win. 

2.2. Hybrid Monte Carlo 

A good deal of progress in integrating the equa- 
tions of motion along an HMC trajectory should 
allow, when combined with the non-hermitian 
solvers of section 1.1, a gain 0(10) in efficiency. 
Here are the main ingredients: 

1) Introduce different "time steps" when comput- 
ing the force coming from the gauge fields and 
that coming from the external boson fields [13]. 
The former can be evaluated cheaply at frequent 
intervals. This will reduce the overall discretiza- 
tion error in the molecular dynamics integration. 
Ref. [14] claims savings by a factor ~ 2. 

2) Use previous solutions of Ax = fo at ear- 
Her time-steps on the trajectory. The simplest, 
widely used extrapolation, consists in choosing as 
a starting guess at step i: Xoi = 2a;,_i — a;,_2- A 
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Figure 6. Logarithm of the change in the gauge 
field caused by a small noise at time 0, versus 
time. From [14]. 

more efficient approach [15] minimizes || Ax — b\\ 
over the span of the earlier solutions, achieving 
another savings ~ 2. My understanding is that 
the system Ax = b could be, in addition, "deflat- 
ed" from the directions already explored, leading 
to faster convergence still. 

3) Use an adaptive time step. This is a very 
promising development, since HMC trajectories 
at low quark mass tend to bounce off energy barri- 
ers caused by a small determinant separating dif- 
ferent topological sectors. It would make sense to 
reduce the step size during those "sharp curves" , 
and accelerate afterwards. Paradoxically, it is 
possible to implement a variable step-size while 
preserving reversibility of the molecular dynam- 
ics evolution [16]. Consider an elementary evolu- 
tion operator = e^(9)^/2e^(p)^e^(«)^/2^ where 
A and B act on conjugate variables p and q, and 
T is the step size, r can be varied adaptively as a 
function of some error e(T, {p, q}{t)) provided 

• the error function is symmetrized, ie. 
e,s(r, {p,q}{t)) = 

e(r, {p, q}it)) + e(-r, Eri{p, 

• at each step the bound on the error must be 
saturated. Thus r is the solution of a nonlinear 
equation es{T, {p,q}{t)) = tolerance. 

These 2 conditions guarantee that the trajectory 
can be retraced step by step, by reversing the 
momenta. This algorithm, tested on the Kepler 
problem with spectacular results, awaits imple- 
mentation for QCD. 

In spite of these promising developments, cau- 
tion is required before jumping to larger lattices 
and smaller quark masses. Roundoff errors be- 
come sizeable on global sums and inner products, 
especially on 32-bit machines. The molecular dy- 



namics integration amplifies these errors expo- 
nentially, as best shown in [14] who measured the 
Lyapunov exponent characteristic of this chaotic 
behavior (sec Fig. 6). HMC relics on long trajec- 
tories, ~ 1 /m, where m is the smallest mass in the 
system, to maintain a dynamical exponent z = 1, 
and on reversibility to converge to the correct dis- 
tribution. Roundoff errors introduce dissipation 
and spoil reversibility. For current trajectories of 
~ 100 steps or more, at the smallest quark masses 
and on the largest lattices, lack of reversibility be- 
comes a serious problem [17]. 

2.3. Kramer's algorithm 

This algorithm, proposed by Horowitz [18] and 
also known as L2MC [19], is a variant of HMC 
based on the second-order Langevin equation: 
q = dH/dp 

p = -dH/dq -^p+ 1/27^ 

where 7 is the viscosity coefficient and t] a 
Gaussian noise. The guiding trajectory of HMC 
now includes some dissipation. The change to 
the usual leapfrog integration scheme is trivial: 
one intercalates an irreversible step p e~'''^p + 
\/l — e~2T'^r?. In addition, all momenta must be 
reversed after rejection by the Metropolis test. 

The original motivation was that the extra tun- 
able parameter 7 could be adjusted to acceler- 
ate the dynamics. The optimal value of 7 has 
been studied for free field in [19]: not surpris- 
ingly, jopt ~ 1/™) so that momenta are effec- 
tively refreshed after a time ~ 1/m as in HMC. 
Neither algorithm appears clearly more efficient 
[14]. Kramer's algorithm, however, maintains a 
dynamical exponent z = 1 for arbitrary short tra- 
jectories. 

I see two advantages in favor of Kramer's algo- 
rithm: 

• Since the dynamics is intrinsicafiy dissipative, 
the additional dissipation due to roundoff errors is 
negfigible. Problems of irreversibility unearthed 
with HMC can safely be ignored. 

• By making trajectories shorter, more configu- 
rations will be generated for the same computer 
cost. They wfil not be completely independent, 
but for some observables with short integrated 
autocorrelation time (like glueballs), they effec- 
tively will be. 
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2.4. Liischer's method 

Liischer's original proposal [20] has now been 
studied thoroughly; it has been modified and im- 
proved a great deal. The two ingredients are: 

• approximate det{A) by det^^ {P{A)), where 
z P{z) « 1 over the whole spectrum of A (see 

Fig.l, with (Jmini^) « mq) 

• factorize P(A) oc Ofcl^ ~ ^k), group the factors 
into positive pairs, and express each as a bosonic 
Gaussian integral. Thus 



{A-Zk)(j>k 



(4) 



The approximation is controlled by the degree 
n of P (ie. the number of bosonic species), and 
converges exponentially: 
\\ P(A) - 1| < ci e-"^ " / m,^ 
with Ci,C2 ~ for all eigenvalues A of 

A. Therefore it is sufficient to increase n as 
m~^LogV with the quark mass niq and the vol- 
ume V of the lattice to keep the quality of the 
approximation constant. 

Early results showed disappointingly slow 
Monte Carlo dynamics [21]. In fact one observes 
that, if the bosonic fields (pu are frozen, the gauge 
fields decorrelate very quickly. Thus the slow 
dynamics are governed by the correlation length 
of the bosonic terms. Its maximum is 0(m^^). 
Calling z the dynamical exponent of the bosonic 
Monte Carlo, one gets for the complexity of the 
algorithm: 

Nb. of bosonic fields oc 
Autocorrelation time oc m^^ 
Work oc m^^^^ 

Although over-relaxation has been used, all re- 
sults are consistent with z = 2. Since the bosonic 
action is a simple Gaussian, there is some hope 
that cluster or multigrid Monte Carlo will provide 
a reduction in z. Even with z = 2, Liischer's algo- 
rithm behaves better asymptotically than Hybrid 
Monte Carlo, where the work grows at least as 

^5/4^-13/4 [22]. 

A large reduction in the work and in the num- 
ber of bosonic fields comes by implementing for 
P{A) the same ideas which have proven success- 
ful for the quark propagator calculation: even- 
odd preconditioning [23], non-hermitian variant 
[24] (which also allows the simulation of an odd 
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Figure 7. Metropolis acceptance versus number of 
bosonic fields in the exact, non-hermitian Liischer 
algorithm; the lattice is 4^; (3 = 5. From [29]. 

number of fiavors) . 

Most significant, however, has been the sugges- 
tion to make the algorithm exact with the adjunc- 
tion of a Metropolis test [25]. The acceptance for 
a transition U ^ U' should then be 

. det-^A{U')P{A(U')) ^ 
' det-^AiU)PiAiU)) ' 

Early efforts aimed at evaluating the determi- 
nant ratio exactly with a Lanczos process [25,27]. 
This approach was plagued by roundoff errors and 
by a cost scaling like V'^. Instead, it is suffi- 
cient to estimate this ratio stochastically [24,26] 
with one Gaussian vector t] as e"**' w-i)ri ^j^-j^ 
W = [A{U')PiA{U'))]-^AiU)PiA{U)). A proof 
of detailed balance is given in [28] . Each Metropo- 
lis test then entails the solution of a linear system; 
one can show that the associated cost, measured 
in update sweeps, remains constant as V and 
increase. 

As an example of the results currently ob- 
tained with this algorithm [29], Fig. 7 shows the 
Metropolis acceptance as the number of fields n 
is varied, for 2 different quark masses. 

A comparison of efficiency with HMC must be 
done carefuhy. Fig. 8 shows some very prehmi- 
nary data, presented on a scale which attempts 
to remove volume efTects. ^^'brk is measured in 
matrix-vector products per independent configu- 
ration; for lack of better data, configurations for 
which the plaquette is decorrelated are considered 
independent. The impression I want to convey 
from this figure is that the exact, non-hermitian 
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Figure 8. The number of matrix-vector products 
necessary to decorrelate the plaquette, divided by 
V^^^,is shown function of the lowest singular 
value of the Wilson matrix, for HMC (X) [17,14] 
and Liischer (+) [29] simulations. 

variant of the Liischer algorithm is a reasonable 
alternative to HMC. 
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